Biexcitons in two-dimensional systems with spatially separated electrons and holes 
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The binding energy and wavefunctions of two-dimensional indirect biexcitons are studied ana- 
lytically and numerically. It is proven that stable biexcitons exist only when the distance between 
electron and hole layers is smaller than a certain critical threshold. Numerical results for the biex- 
citon binding energies are obtained using the stochastic variational method and compared with the 
analytical asymptotics. The threshold interlayer separation and its uncertainty are estimated. The 
results are compared with those obtained by other techniques, in particular, the diffusion Monte- 
Carlo method and the Born-Oppenheimer approximation. 



I. THE PROBLEM AND MAIN RESULTS 

The physics of cold excitons — bound states of elec- 
trons and holes in semiconductors — has attracted much 
attention recently. Cooling the excitons has become pos- 
sible by confining electrons and holes in separate two- 
dimensional (2D) quantum wells, which greatly increases 
their lifetime. A number of intriguing phenomena has 
been demonstrated for such "indirect" excitons, includ- 
ing long-range transport, 1-6 macroscopic spatial order- 
ing, 4 and spontaneous coherence. 7 Theoretical work on 
these phenomena is ongoing, see Ref. 8 for review. Fur- 
ther progress in this field requires an improved under- 
standing of exciton interactions. 

Despite being charge neutral, indirect excitons possess 
a dipole moment ed, where d is the separation of the 
electron and hole quantum wells. As a result, interaction 
of two excitons at large distances r is dominated by their 
dipolar repulsion, 



V(r) = 



(1) 



where k is the dielectric constant of the semiconductor. 
At short distances exchange and correlation effects are 
also important. The interaction may even become at- 
tractive over a range of r. In this case two excitons can 
form a bound state — a biexciton. The corresponding 
binding energy is defined by 



Eb = 2Ex — E 



xx • 



(2) 



where Ex and Exx are the ground-state energies of the 
exciton and biexciton, respectively. 

While observations of biexcitons in single quantum 
well structures (d = 0) have been described multiple 
times, 9-16 no such reports exist for the d > case. A 
recent theoretical work 17 has attributed the lack of ex- 
perimental signatures of indirect biexcitons to extreme 
smallncss of their binding energies. In this paper we ver- 
ify and improve all previously known estimates of Eb- 
In particular, we show that Es{d) is positive, i.e., the 
biexciton is stable, only for d smaller than some critical 
value d c , see Fig. 1. Typical experimental parameters 8 ' 18 
fall on the d > d c part of the diagram. 



In our calculations we adopt the simplifying assump- 
tion that the effective masses m e and rrih > m e of elec- 
trons and holes are constant and isotropic. We also treat 
the quantum wells as 2D layers of zero thickness. We 
find it convenient to measure distances in units of the 
effective electron Bohr radius and energies in units of the 
effective Rydberg, 



p 1 6 

2 K<2„ 



(3) 



respectively. With these conventions, the four-particle 
system of two electrons and two holes is described by the 
Hamiltonian Hxx — T + U, where 



U = 



T = Ti- 
2 

|ri - r 2 



T 2 



Tj = -V 
2 



(4) 



|Ri- 
v(r,d) 



R 2 



-]>>(!-, -R,,d), (5) 



(6) 



V\t\ 2 + <p 

Here and Ri are 2D coordinates of the electrons and 
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FIG. 1: Critical interlayer separation vs. the electron-hole 
mass ratio. The circles are our results. The squares are from 
Ref. 24. The triangles correspond to d above which J5s(d) 



drops below 10 ; 
imental practice. 



Ry e , making biexcitons irrelevant in exper- 
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the holes, respectively, Vj = d/drj, and 
a = m e /m h 



( 7 ) 



is the mass ratio. Similarly, the single-exciton Hamilto- 
nian is 



H x =T 1 +v{r 1 -Ri,d). 



(8) 



The problem is characterized by two dimensionless pa- 
rameters: d and a. The case of d = (direct excitons) 
has been studied extensively. 19-21 In contrast, high accu- 
racy calculations of Eb for d > have been carried out 
only in the aforementioned Ref. 17. The authors of that 
work have employed the diffusion quantum Monte-Carlo 
method (DMC). Away from d = 0, they were able to fit 
their results for a = 1 and a = 1/2 to the exponential: 



Enid) ~ ae 



-0d 



(9) 



This result is surprising. Equation (9) seems to imply 
that the biexcitons are stable at any d, i.e., d c = oo. On 
the other hand, physical intuition and previous approxi- 
mate calculations 22,23 suggest that d c should be finite. A 
more recent work 24 has reached the same conclusion. In 
this paper we present rigorous analytical arguments and 
essentially exact numerical results proving that d c < 1 at 
all <t, see Fig. 1. (Due to electron- hole symmetry, it is 
sufficient to consider < a < 1.) 

Since d c is finite, the interpolation formula (9) must 
overestimate the binding energy at large d. We show 
that near the biexciton dissociation threshold, 



d< D. 



(10) 



where D ~ 1 for a ~ 1 and D ~ exp(— a 1 / 2 ) for a <C 1, 
function Eg(d) behaves as 



£ e 



-D/(d c -d) 



(11) 



This equation resembles the well-known expression for 
the energy £ of a bound state in a weak 2D potential 
V(r). Such a state exists if 



W 



M 
2wh 2 



d 2 rV{r) < 0. 



(12) 



where M is the mass of the particle. Near the threshold 
W — > one finds 25 



|e| oce- 1/|w| , |W|««1 



(13) 



The exciton-exciton interaction potential V(r) in general 
does not satisfy the condition of the perturbation theory 
V(r)r 2 < h 2 /M, withM = m e +m h . Therefore, Eq. (13) 
does not literally apply here. Nevertheless, the physical 
origins of the exponential dependence in Eqs. (11) and 
(13) are the same, see Sec. IIB. 

We verify and complement the above analytical re- 
sults numerically using the stochastic variational method 
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FIG. 2: Binding energy vs. the distance between the quantum 
wells for the mass ratios a = 1 and 0.5. Our results are shown 
by the solid lines. The dots are from the Ref. 17. 



(SVM). 26 The SVM has proven to be a powerful tech- 
nique for computing the energies of few-particle sys- 
tems. 27 For example, it has given the best estimates of 
Eb for direct biexcitons, 19,20 d = 0. Our calculations are 
largely in excellent agreement with those of Ref. 17, see 
Fig. 2 and Table I. Thus, Eq. (9) is certainly useful as 
an interpolation formula for not too large d. However, 
near the estimated d c , our results favor Eq. (11) over 
Eq. (9). Since the SVM is variational, we can be sure 
that it is more reliable when it gives a larger Eb than 
other methods. 

The remainder of the paper is organized as follows. 
In Sec. II we derive a few analytical bounds on Eb and 
the asymptotic formula (11). Numerical calculations are 
presented in Sec. III. Section IV is devoted to discussion 
and comparison with results in previous literature. Some 
details of the derivation are given in the Appendix. 



II. ANALYTICAL RESULTS 

In this section we approach the biexciton problem by 
analytical methods. Since the exact solution seems out 
of reach, the best one can do is to consider certain limits 
where suitable control parameters exist. Below we ex- 
amine three of them. First, we study large-ti excitons. 
We prove that they cannot bind into a stable biexciton. 
Second, we consider the immediate vicinity d c — d -C 1 
of the dissociation threshold d c . We derive the asymp- 
totical formula for the binding energy, Eq. (11), which is 
valid for arbitrary a. Finally, we analyze the case a <C 1. 



A. Exciton interaction at large d 

The absence of stable biexcitons at large d is due to 
the lack of binding in the classical limit, which is realized 
at such d. Indeed, if we temporarily change the length 
units to d and energy units to e 2 / nd, then the potential 



3 



energy U in Eq. (5) becomes d-independent while the 
kinetic energy T acquires the extra factor a e /d <C 1 com- 
pared to Eq. (4). Hence, the potential energy dominates. 
A rigorous proof that d c < oo can be constructed by 
dealing with the quantum and many-body aspects of the 
problem separately. The many-body part is handled at 
the classical level. Thereafter the quantum corrections 
are included. With further analysis, both parts of the 
argument can be reduced to simpler problems for which 
controlled approximations exist. 

Since the Earnshaw theorem does not apply in 2D, the 
absence of a stable classical biexciton is not immediately 
obvious. However, we verified it following these steps. 
The classical ground-state is the global minimum of the 
potential energy. We can do the minimization over the 
electron positions ri and r 2 first. Let R be the distance 
between the holes, 

R = Ri R 2 , (14) 




1 2 3 4 5 

R/d 



FIG. 3: Main panel: ground-state energy U m i n vs. the sepa- 
ration R of holes for a pair of classical excitons. In this state 
all four charges are on the same straight line. Inset: in-plane 
distance between nearest electrons and holes vs. R. 



then the energy function to minimize is (in the original 



TABLE I: Biexciton binding energies in units of Ry e from the 
previous ("DMC", Ref. 17) and present ("SVM") work. 



d(tte) 


a — 1 
DMC 


SVM 


d(tte) 


cr = 0.5 

DMC 


SVM 


0.00 


0.3789 


0.3858 


0.00 


0.5381 


0.5526 


0.02 


0.3084 


0.3089 


0.01 


0.4443 


0.4450 


0.04 


0.2538 


0.2546 


0.03 


0.3695 


0.3689 


0.06 


0.2118 


0.2133 


0.04 


0.3104 


0.3109 


0.08 


0.1794 


0.1807 


0.06 


0.2639 


0.2649 


0.10 


0.1532 


0.1542 


0.07 


0.2265 


0.2275 


0.12 


0.1315 


0.1324 


0.09 


0.1956 


0.1966 


0.14 


0.1135 


0.1141 


0.11 


0.1696 


0.1707 


0.16 


0.0982 


0.0986 


0.12 


0.1477 


0.1487 


0.18 


0.0851 


0.0855 


0.14 


0.1291 


0.1299 


0.20 


0.0738 


0.0742 


0.15 


0.1130 


0.1136 


0.22 


0.0640 


0.0644 


0.17 


0.0989 


0.0995 


0.24 


0.0556 


0.0559 


0.18 


0.0865 


0.0872 


0.26 


0.0483 


0.0485 


0.20 


0.0757 


0.0764 


0.28 


0.0418 


0.0420 


0.21 


0.0663 


0.0670 


0.30 


0.0361 


0.0363 


0.23 


0.0580 


0.0586 


0.32 


0.0311 


0.0313 


0.24 


0.0507 


0.0512 


0.34 


0.0267 


0.0270 


0.26 


0.0443 


0.0447 


0.36 


0.0229 


0.0231 


0.27 


0.0385 


0.0389 


0.38 


0.0195 


0.0197 


0.28 


0.0333 


0.0337 


0.40 


0.0165 


0.0167 


0.30 


0.0286 


0.0291 


0.42 


0.0140 


0.0141 


0.32 


0.0241 


0.0250 


0.44 


0.0117 


0.0118 


0.33 


0.0200 


0.0214 


0.46 


0.0096 


0.0097 


0.34 


0.0165 


0.0182 


0.48 


0.0078 


0.0079 


0.36 


0.0135 


0.0154 


0.50 


0.0063 


0.0065 


0.38 


0.0112 


0.0129 


0.52 


0.0051 


0.0052 


0.39 


0.0096 


0.0107 


0.54 


0.0040 


0.0040 


0.41 


0.0087 


0.0087 


0.56 


0.0030 


0.0031 


0.42 


0.0076 


0.0071 


0.58 


0.0021 


0.0023 


0.44 


0.0064 


0.0056 


0.60 


0.0013 


0.0017 


0.45 


0.0051 


0.0044 


0.62 


0.0007 


0.0012 


0.47 


0.0039 


0.0033 


0.64 


0.0002 


0.0007 


0.48 


0.0027 


0.0024 



units convention) 

1 1 .7=1,2 

t=±R/2 

It can be shown that for all R the lowest energy 
is achieved when the in-plane coordinates of the four 
charges fall on a straight line, see Fig. 3. Forming a cross 
is the only other viable alternative, but it always has a 
higher energy. For the linear geometry of the system, nu- 
merically exact results for U m i n (R, d) = min riiI . 2 Ur are 
obtained trivially. The plot of V C \(R) = U m i n (R,d) + 
(4/d) is shown in Fig. 3. This combination can be 
thought of as the classical limit of the exciton interaction 
potential V(R). Function V c \ monotonously decreases 
with R and achieve its global minimum at R = oo. This 
means that classical excitons do not form a bound state. 

At large R, function V c i(R) follows the dipolar inter- 
action law (1) with the quadrupolar, etc., corrections: 

Id 2 3d 4 / d 6 \ 

v ^ R ^ = w-2i +0 {w)' Ryyd - (16) 

Quantum corrections due to the zero-point motion about 
the classical ground state are not able to compete with 
the dipolar repulsion when d is large, see Appendix A. 
Therefore, there is a critical d c = d c {a) above which a 
stable biexciton does not exist. 



B. Binding energy near d c 

In this subsection we examine the biexciton state near 
the dissociation threshold d c for arbitrary a. It is easy 
to understand that in this regime the biexciton orbital 
wavefunction 'F should have a long tail extending to large 
distances away from the center of mass of the the sys- 
tem. Inside of this tail the configurations of electrons 
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and holes resemble a pair of well-separated individual 
excitons. Therefore, at r > 1, where r is the distance 
between the centers of mass of two such excitons, VP takes 
the asymptotic form 

* = [i + (-iyp 12 ] *(r) n M*j R,-) , (17) 

3=1,2 



(18) 



Here s is the total electron spin, <j) a is the ground-state 
wavefunction of a single exciton with mass ratio a, and 
operator P i2 exchanges ri and r 2 . Let us assume, for 
simplicity, that holes are spin- 1/2 particles. Then the 
wavefunction $ of the relative motion must have the par- 
ity r) = (-l) s+s <l>(r), where S is the total spin of 
the holes. Our goal in this subsection is to determine the 
behavior of <J> at large r and use it to derive Eq. (11). 

We proceed, as usual, by expanding $ into partial 
waves of angular momenta m (m and s + S must be 
simultaneously odd or even). The equation for the radial 
wavefunction Xm{ r ) reads 



ld_ dxm 
r dr dr 



x 2 + (j,V(r) + -5- 



Xm=0, (19) 



where x and /i are defined by 



1 + q 

2(7 



(20) 



At small distances, potential V(r) is either ill-defined 
or complicated, but for r 3> d it obeys the dipolar law 
V(r) = 2d 2 /r 3 [Eq. (1)]. From this, it is easy to see that 
/iV(r)r 2 < 1 at r > 6 with b given by 



b = 8fid 2 . 



(21) 



At such r the potential energy V acts as a small perturba- 
tion. 25 Therefore, Xm(r) coincides with the wavefunction 
of a free particle, 



Xm(r) = c\K m (xr) , r > b . 



(22) 



Note that b is either of the order or much larger than d 
because fi > 2 and d ~ <i c ~ 1. 

Sufficiently close to the critical d, we have x <C 1/6. 
In this case there exists an interval of distances 6 <C r <C 
& 1 / 3 x -2 / 3 where we can drop the term x 2 in Eq. (19) 
compared to fiV(r). After this, Eq. (19) admits the so- 
lution 



(23) 



where hm{z) and K<2 m (z) are the modified Bessel func- 
tion of the first and the second kind, respectively. 28 The 
unit coefficient for / 2m (z) and the factor of (—4) in front 
of c 2 are chosen for the sake of convenience. The ground- 
state solution is obtained for m — 0. Using the asymp- 
totic expansion 28 of I and K in Eqs. (22) and (23), and 




demanding them to be consistent with one another, we 
find for m = and b <C r <C x~ 1 : 



Xo = 1 - 2c 2 



C2 



-27 



+ 
1 



67 + 2 ln(6x/8) ln(E /E B ) 



where 



En — 



e 6 T V 1 + cr 



1 



(24) 
(25) 

(26) 



Here 7 = 0.577 ... is the Euler-Mascheroni constant. 28 
Equation (24) specifies the boundary condition to which 
the solution for xa m the near field, r < b, must be 
matched. 

At d = d c , both x and c 2 vanish. Wavefunction Xo( r ) 
at small x can be viewed as the wavefunction for d = d c 
perturbed by the small change in the boundary condition 
in the far field, r >b, and by another perturbation, 

** + l*V(r)\i e , 

in the near field, r < b. To the first order in these per- 
turbations we have 



Er 



-Ac2 + B{d,x 2 ), 



(27) 



where A is a constant and B is a smooth function subject 
to the condition B(d c ,0) = 0. Expanding B to the first 
order in d c — d and x 2 , we arrive at the transcendental 
equation for Eb'- 



dB\ A 

1 -^ 2 ) Eb + Me /e b ) 



= -^(d c -d). (28) 



The solution cannot be written in terms of elemen- 
tary functions. However, the logarithmic term gives the 
sharpest dependence on Eg. Hence, at small Eg the first 
term on the left-hand side of Eq. (28) can be dropped. 
Now this equation can be easily solved to recover Eq. (11) 
with 



C ~ dd 



(29) 



The coefficients A and C must be determined from the 
solution of the inner problem. For a < 1 part of this 
task can be accomplished analytically, as explained later 
in this section. For j ~ 1 a numerical solution, such 
as the one discussed in Sec. Ill, seems to be the only 
alternative. 

Our results comply with a general theorem, 29 which 
states that in the asymptotic limit k = ix — > the scat- 
tering phase shift S(k) satisfies the equation 



(n/2) cot S(k) = ln(fc/2) + f(k 2 ) , 



(30) 



where f(z) is some analytic function. This theorem is 
valid for a general short-range potential in 2D. For a 
bound state cotS(k) should be replaced by i, leading to 



lnk/A/2 +/Mb) = 0, 



(31) 
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which is in agreement with our Eq. (28). Our derivation 
has the advantage of showing that the proper dimen- 
sionless combination in the argument of the logarithm is 
s/Eb/Eo and that the asymptotic behavior (11) is real- 
ized at Eb <C E n . 



C. Binding energy for small mass ratios 

Although the electron-hole mass ratio is not truly small 
in typical semiconductors, it is interesting to examine the 
case a < 1 from the theoretical point of view. At such 
a the exciton interaction potential V can be meaning- 
fully defined at all distances using the Born-Oppenheimer 
approximation (BOA). 30,31 In addition, the radial wave- 
function can be computed everywhere with accuracy 
O(a). 

The distance r between excitons is no longer a phys- 
ically reasonable variable when the four particles ap- 
proach each other closely and their partitioning into ex- 
citons becomes ambiguous. In the BOA this problem 
is mitigated by selecting R — the distance between the 
heavy charges — to be the radial coordinate of choice. 
The ground-state biexciton wavefunction is taken to be 



tf = xCRMR,ri,r 2 ), 



(32) 



where ip is the ground-state of two interacting electrons 
subject to the potential of two holes fixed at positions 
Ri, 2 - ±R/2: 

Hboav = [-V? - V 2 . + U R {n,r 2 )] <p = U BO Af ■ (33) 

Here Uboa(R) is the corresponding energy. In turn, x(i?) 
is found from 

-^f^ R m +fx[UBOA{R) ~ EBOA]x=0 - (34) 

The BOA is known to have 0{a) accuracy. In principle, it 
can be systematically improved. 32 However, since below 
we will be solving Eq. (34) by means of the quasiclassical 
approximation, which itself is known to be accurate only 
up to 0(c), this is unwarranted. 

Dropping all inessential 0{a) terms, we can simplify 
Eq. (34) as follows: 



1 Ao* 
R dR dR 



+ [^ + fiV(R)] X = 0, (35) 
V(R) ee Uboa(R) - CW(oo) • (36) 



Our task is to solve this equation with boundary condi- 
tions |x(0)| < oo at the origin and 



x(R) 




4c 2 K 




(37) 



at b < R < & 1/3 xr~ 2/3 , with c 2 given by Eq. (25). 

We reason as follows: in order to have a bound state, 
potential V(R) must be negative over some range of R. It 



can be shown that this occurs in a single contiguous inter- 
val, see Fig. 4 and Sec. III. Inside of this interval there is a 
classically allowed region, /iV(R) < — x 2 , where function 
x(i?) reaches a maximum. As we approach the dissocia- 
tion threshold, this region shrinks. Near the threshold it 
becomes very narrow, so that the quadratic approxima- 
tion 



(iV(R) 



\»V"{R 



R~){R-R^ 



(38) 



becomes legitimate. Here i?_ and R + are the turning 
points. To construct the desired solution we simply need 
to match x(i?) in the classical region, i?_ < R < R + , 
inside the tunneling region, R+ and in the far 

field, R 3> b. Details of this calculation are outlined in 
Appendix B. The result is 



A = 4(%F") 1/2 cxp(-2S ), 
dRy/V{R) , 



Sq 



(39) 
(40) 
(41) 



B = \V(R )\-V^V 77 , 

where Rq = (R+ + i?_)/2 is the point where V(R) has 
the minimum. 

Equations (29) and (39) imply that the coefficient D 
in Eq. (11), and so the range (10) of d where Eq. (11) 
applies are proportional to the exponentially small factor 
e~ 2S ° at a <C 1. We expect that D grows with a and by 
extrapolation, reaches a number of the order of unity at 
<r~ 1. 

A few other properties of function d c (a) can also be 
deduced analyticaly. For example, Eq. (41) implies that 



^c(O) — d c (cr) oc \fa , a <C 1 



(42) 



Hence, d c (a) has an infinite derivative at a — and so 
initially decreases with a. At some cr, however, d c (a) 
must start to increase. Indeed, due to the electron-hole 



\V(R) 










\ x(R) 




















V-hH — Y-f- 

\ R / 


— i — 

d 


1 

b 


x- 1 


R 





FIG. 4: Sketch of the interaction potential V(R) and the 
exciton wavefunction x(R) f° r the Born-Oppenheimer limit 
a < 1. 
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symmetry, the combination d c (a)/(l + a) must have a 
vanishing derivative 33 at a = 1. Therefore, 

d' c (l) =d e (l)/2 > 0. (43) 

Finally, we have a strict upper bound 33 

4(a) < (1 + ^)4(0). (44) 

All of these properties are borne out by our Fig. 1. Still, 
a purely analytical solution of the biexciton problem does 
not appear to be possible at any a. In the next section, 
we approach it by numerical calculations. 

III. NUMERICAL SIMULATIONS 

In order to verify our analytical predictions and other 
results in the literature 17,24 , we have carried out a series 
of numerical calculations using the SVM. To implement 
this method we customized the published SVM computer 
code 34 for the problem at hand. In the SVM one adopts a 
nonorthogonal basis of correlated Gaussians in the form 27 

G n = exp (-~iJB n x\ , (45) 

from which a variational wavefunction of given electron 
and hole spins (S and s, respectively) is constructed: 

* = A[G n ({T v })? s , a ]. (46) 

Here x is a 3 x 1 vector of Jacobi coordinates (linear com- 
binations of differences in particle coordinates in which 
the kinetic energy separates), B„ is a positive-definite 
3x3 matrix, A is the antisymmetrizer, and Ts, s is the 
spin wavefunction. All our SVM calculations are done 
for the spin-singlet state S = s = 0. Note that G n corre- 
sponds to the zero total momentum of the system. 

The number of basis states is grown incrementally un- 
til the energy is converged or the prescribed basis di- 
mension (typically 700) is reached. At each step a new 
quadratic form B„ is generated randomly. If adding the 
corresponding function G n to the basis improves the vari- 
ational energy significantly, this G n is kept; otherwise, a 
new B„ is generated by varying some of its matrix ele- 
ments. Details can be found in Refs. 27 and 34. 

Our numerical results for a — 0.5 and a = 1 are given 
in Table I and plotted in Fig. 2. In Fig. 5 we replot the 
binding energy Eb for a = 1 in a form suitable for testing 
Eq. (11): 

1 = d c -d (d c - df 
ln{E /E B ) D D x ' 1 ' 

Here we take into account one more term in the Taylor 
expansion of the right-hand side of Eq. (27) compared 
to Eq. (28). Extrapolation of the data to Eb = gives 
us d c . The uncertainties in this parameter are estimated 
by imposing a 95% confidence level on the fit coefficients 




d{a e ) 

FIG. 5: Logarithmic plot of the biexciton binding energy as 
a function of d for a = 1. Our results are shown by the filled 
symbols; the open circles are from Ref. 17. The thicker line is 
the fit to Eq. (47), which yields d c = 0.87 ± 0.01 with a 95% 
confidence level. The other line is Eq. (9) with a and /3 from 
Ref. 17. 

d Cl D, and D\. The same procedure has been applied 
to several other mass ratios in the interval 0.1 < a < 1. 
The results for d c are shown in Fig. 1. Their comparison 
with other results in the literature will be addressed in 
Sec. IV. 

At a < 0.1 the range (10) of d where Eq. (11) ap- 
plies is exponentially small. Even with our highly accu- 
rate numerical method we were not able to probe this 
range. Thus, we assumed that the nonanalytical correc- 
tion Ac2(E B ) is undetectable on the background of Eb 
in Eq. (27), so that our numerical results for Esid) at 
such a are dominated by the regular contribution 

E B = C(d c -d) + C 1 (d c -d) 2 + ... (48) 

Accordingly, at a < 0.1 we deduced d c from the fit 
of Es(d) to a quadratic polynomial. Additionally, we 
confirmed that at a = 0.2 the two fitting procedures 
give similar results: d c = 0.59 ± 0.01 per Eq. (47) vs. 
d c = 0.58 ± 0.01 per Eq. (48). 

Finally, we have computed the electron and hole den- 
sities in the biexciton as a function of their distance from 
the center of mass. Examples are presented in Fig. 6 for 
d = 0.0 and d = 0.3. In the latter case the particles are 
on average further away from the center of mass. The 
same trend is also seen in the average root-mean square 
separations between various particles, which are plotted 
in Fig. 7. Their accelerated growth with d occurs be- 
cause the biexciton becomes less bound and eventually 
dissociates. 



IV. DISCUSSION 

Let us compare our results with previous theoretical 
work. Early studies of the biexcitons based on Hartree- 



7 



Fock 23 or Heitler-London 35 approximations provided ini- 
tial evidence for the existence of a finite threshold d c for 
the biexciton dissociation. However, they gave a con- 
siderably lower d c that we find here because these ap- 
proximations did not account for all correlation effects 
essential to the biexciton stability. 

Comparing with more recent calculation 17 of the biex- 
citon binding energies by the DMC technique, we find an 
overall excellent agreement. Still, our SVM occasionally 
slightly outperforms the DMC, see Table I. Furthermore, 
in the SVM the estimate of the ground-state energy de- 
creases at each step, so that the statistical noise is never 
an issue, unlike in the Monte-Carlo methods. Neither the 
SVM nor the DMC is able to compute arbitrarily small 
binding energies; therefore, in order to determine d c , an 
extrapolation to Eg = is necessary. The clarification of 
what extrapolation formula should be used for this pur- 
pose is an important finding of this work. Equation (47) 
represents the true asymptotic behavior in the limit of 
small Eb and indeed describes our numerical results at 
such Eb better than the interpolation formula (9) plotted 
alongside for reference. 

Another recent theoretical work on biexcitons used a 




FIG. 6: (a) Electron and hole density vs. the distance to the 
center of mass in a biexciton with a — 0.5 and d = 0.0. (b) 
Same for a — 0.5 and d = 0.3. 
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FIG. 7: Root mean square of the pairwise distances between 
the biexciton constituents vs. d for a = 0.5 and a = 1. 



Born-Oppenheimer-like approximation to compute the 
threshold interlayer separation d c (a). At all a shown in 
Fig. 1, our d c are higher than those reported in Ref. 24. 
The deviation is much larger than the uncertainty of d c 
from our extrapolation procedure. This is surprising be- 
cause the adiabatic BOA 36 is known to give a strict lower 
bound on the ground state energy 37 while variational 
methods, such as our SVM, give a strict upper bound. 
Since the energy of a single exciton Ex is usually com- 
puted extremely accurately, our binding energies should 
be smaller than those of Ref. 24. Accordingly, their esti- 
mates of d c should exceed ours, not the other way around. 
We suspect that the problem may again be related to 
the manner in which the Eb — > extrapolation was per- 
formed in Ref. 24. In any case, a significant difference 
is seen only at a ~ 1. At small mass ratios, where the 
approximation of Ref. 24 becomes accurate to the order 
0{a), our results are in better agreement. 

Turning to the experimental implications of our theory, 
observations of biexcitons in single quantum well systems 
have been reported by many experimental groups. 9-16 In 
contrast, no biexciton signatures have ever been detected 
in electron- hole bilayers. Let us discuss how this can be 
understood based on our results. 

The first point to keep in mind is that the biexci- 
ton dissociation threshold d c plotted in Fig. 1 is a zero- 
temperature quantity. For the biexcitons to be observ- 
able at finite temperatures, Eb must exceed kT by some 
numerical factor. (As usual in dissociation reactions, 38 
this factor is larger the smaller the exciton density is.) 
The coldest temperatures demonstrated for the excitons 
in quantum wells are T ~ 0.1 K [Ref. 39]. The maximum 
separation between the 2D electron and hole layers 
at which biexcitons are still physically relevant in such 
structures can be roughly estimated from 



E B (d*) = HT 3 Ry e . 

Function d* (a) is plotted by triangles in Fig 
GaAs quantum wells we have 40 a m 0.5, a e - 



(49) 

1. In 
10 nm, 
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and so w 4.5 nm. In comparison, the smallest 
center-to-center separation that has been achieved in 
GaAs/AlGaAs and InGaAs/GaAs quantum wells with- 
out compromising the sample quality is at least twice 
as large. 8 Cold gases of indirect excitons have also been 
demonstrated in AlAs/GaAs structures, 18 in which d is 
smaller, d — 3.5 nm. But the electron Bohr radius is also 
smaller, a e w 3nm, so, unfortunately, the dimensionless 
d is about the same. 

A more serious obstacle to the creation and observation 
of biexcitons is disorder. A rough measure of disorder 
strength is given by the linewidth of the exciton optical 
emissions, which is currently ~ 1 meV, i.e., of the order 
of 0.1 Ry e in GaAs. Eb becomes smaller than this energy 
scale as soon as d exceeds the thickness of a few atomic 
monolayers, see Fig. 2. Actually, if the disorder were due 
to a long-range random potential, it might still be pos- 
sible to circumvent its influence on the measured optical 
linewidth by interferometric methods such as quantum 
beats. 11 ' 13 In reality, a short-range random potential is 
probably quite significant. 

One potentially promising system for the study of the 
biexciton stability diagram is a single wide quantum well 
subject to an external transverse electric field. 41 If the 
well is symmetric and the applied field is zero, we have 
d = 0. A finite field can pull electrons and holes apart, 
leading to d > 0. Of course, for such a structure one 
should recalculate the stability diagram of Fig. 1 by tak- 
ing into account the motion of particles in all three di- 
mensions. 

Although it is challenging to observe the binding of 
free indirect excitons, in experiments they can be loaded 
and held together in artifical traps. 42 We anticipate that 
the SVM can be a powerful tool to study systems of a 
few trapped excitons theoretically, complementing recent 
Monte-Carlo work. 43 

In conclusion, we have obtained the most accurate esti- 
mates to date of the binding energies of two-dimensional 
biexcitons. Future work may include a refined study of 
exciton-exciton scattering 24 and interacting excitons in 
traps. 

This work is supported by the NSF grant DMR- 
0706654. We are grateful to R. Needs for providing us 
with the numerical DMC results from Ref. 17 listed in 
Table I. We thank L. Butov and L. Sham for valuable 
discussions. 



APPENDIX A: RIGOROUS BOUNDS FOR THE 
BIEXCITON BINDING ENERGY 

In this appendix we give a few strict upper bounds on 
Eb, which enable us to prove the nonexistence of sta- 
ble biexcitons at sufficiently large d. The basic logic of 
the proof was outlined in Sec. IIA. Here we provide the 
technical details. 



Our starting bound is 

E B < max E R , (Al) 

R 

where 

Eb = inf spec H^ — inf spec Hr (A2) 

is the binding energy of the two-electron Hamiltonian 
Hr = Tr + Ur whose kinetic term is 

T fl = -(l + a)(V 2 + Vi), (A3) 

and the potential term Ur is given by Eq. (15). The 
Hamiltonian Hr is similar to that of the original problem 
[Eqs. (4)-(17)] except the holes are replaced by static 
charges separated by a given distance R and the electron 
mass is made equal to the reduced electron-hole mass. 

To derive the inequality (Al) we take advantage of the 
well-known theorem that the ground-state energy as a 
concave function in the strength of an arbitrary linear 
perturbation. (This theorem follows from the variational 
principle.) For our purposes we choose the perturbation 
in the form 

"■>-*-{■&;)'■ < A4) 

We add it to the kinetic energy terms with the coefficient 
-<j < t < 1, yielding Tj — ► Tj + tATj. Hamiltonians H 
and Hr are obtained by setting r = and r = — er, 
respectively. 

The perturbation leaves the reduced electron-hole mass 
invariant. Therefore, it does not effect the ground-state 
energy Ex of a single exciton. The energy E xx (t) does 
vary with t and the aforementioned concavity property 
dictates 

Exx(t) > \-^E xx {-a) + p-^E xx (l) . (A5) 
1 + c 1 + a 

Since E xx {— a) = E xx (l) by electron-hole symmetry, 
the right-hand side is equal to E xx (— a) for all r. Con- 
sequently, t = —a gives the largest binding energy and 
we arrive at the inequality (Al). 

If the kinetic energy Tr is discarded, Er becomes equal 
to —V c \(R,d) < 0. We want to ascertain that quantum 
corrections do not change the sign of Er. 

The quantum corrections appear in both E x and E xx . 
The former are well understood. 22 The internal dynamics 
of the exciton in the large-rf case is analogous to that of 
a 2D harmonic oscillator with the amplitude of the zero- 
point motion given by 

(| ri -Ri| 2 > = l 2 , J = d 3 / 4 (l + a) 1 / 4 «d. (A6) 

The corresponding energy correction is 

*-§-^-°(*)- <A7) 
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This result immediately restricts the range of R where 
the stable biexciton may in principle exist. By positivity 
of the kinetic energy, Er < 2Ex — ?7 m i n (i?, d), where 
U m [ n is defined in Sec. IIA. Therefore, Er > may occur 
only at R that satisfy 



Vci(R) > 2E X + 



d 



(A8) 



In view of Eqs. (16) and (A7), R must necessarily be 
much larger than d. 

Choose an arbitrary di such that d <C d\ <C R. By 
definition of C/ m i n , 

U R > [/ mi „(iZ,di) + VV(n) + Vy(r 2 ) , (A9) 

Vy(r) = J2 Kr-t,di)-«(r-t,d)]. (A10) 

t=±R/2 

Accordingly, E R < 2Ex — t^min {R, d\ ) - 2E Y , where Ey 
is the ground-state energy of a single electron subject to 
the potential Vy (r) of four out-of-plane charges. This po- 
tential has the shape of two symmetric wells separated by 
the distance R. The amplitude of the zero-point motion 
in each well is again I R. Therefore, the energy shift 
due to tunneling between the wells is exponentially small. 
(A rigorous upper bound can be given. 44 ) Furthermore, 
potential Vy near the bottom of each well coincides with 
that of a single exciton up to a constant 



AVy = Vy 



R 



_ 2 _ 2_ d\-d 2 
~d~d~ 1 + i? 3 



Hence, E Y = E x + AVy and 



E < 2d " 



Vcl (R, dl )- 2 -§ 



(All) 



(A12) 



In these formulas we have dropped subleading terms 
o(l 2 /d\), o(df/R 5 ), etc. With the same accuracy the 
bracket in Eq. (A12) vanishes [cf. Eq. (16)], so that we 
arrive at the result Er ~ —V C \{R, d). This simply means 
that at large d all quantum corrections to Er are para- 
metrically smaller than the direct dipolar repulsion of 
the two excitons. Therefore, Er < at all R, so that 
Eb < 0, and the proof is complete. 



APPENDIX B: RADIAL WAVEFUNCTION FOR 
SMALL MASS RATIOS 

In this appendix we show how the suitable solution 
of Eq. (35) can be constructed within the quasiclassical 
approximation. The necessary connection formulas are 
derived by asymptotic matching with two exact solutions 
at small and large R. 

It is convenient to define the rcscalcd wavefunction 
ip(R) = x(-R) VR- From Eq. (35) we find that tp satisfies 
the equation 



r - x 2 + ^v(r) 



0. 



(Bl) 



This equation has two linearly independent quasiclassical 
solutions 



V±(i?) = 



1 



exp ( ± [S(R) - 5(6)] ) , (B2) 



VOW) 
where Q and S are given by 

R 

Q(R) - V* 2 + t*V{R) , S(R) = J dpQ(p). (B3) 

R+ 

The subtraction of the ^-independent term S(b) in the 
exponentials amounts to multiplying ip± by unimportant 
constants. This is done purely for the sake of conve- 
nience. The reason for omitting the l/4i? 2 term in the 
formula for Q is more subtle. It is explained in detail in 
Rcf. 45. 

In the following we assume that x <C 1/6, in which case 
there exists a broad interval d <C R <C 6 where potential 
V{R) is dominated by the dipolar repulsion (1). In this 
interval, /jV(R) ~ 6/4i? 3 > k 2 \ therefore, 



i, ± (R) 



b R 



1/4 



exp 




(B4) 



Using the asymptotic expansion formulas 28 for I and 
Kq, it is easy to see that the following linear combination 



^(R)^—4,_(R) 



2y^ 



C2^+{R) 



(B5) 



of the quasiclassical wavefunctions (B4) smoothly 
matches with the exact solution (37) at d <C R <C 6. 
This is our first connection formula. It is crucial for this 
derivation because in the intermediate range of distances 
6 <C R <C x the quasiclassical approximation breaks 
down. (It is invalidated by the sharp decrease of V with 
R.) In that region x{R) — ip/VR exhibits a slow log- 
arithmic falloff (24) instead of the algebraic decay sug- 
gested by Eq. (B4). As explained in Sec. II, the nonan- 
alytical behavior (11) of the binding energy is precisely 
due to this logarithmic falloff. 

To finish the calculation we need a second connection 
formula between \ given by Eq. (B5) and the same func- 
tion near the classical turning point R+. To find it we 
take advantage of the exact solution for the harmonic os- 
cillator potential (38) in terms of the parabolic cylinder 
function 28 



tjj cx D e _ 1/2 (-V / 2a:) . 



R — Rt 



o 



(B6) 



Here Ro = (i?+ + i?_ ) /2 is the point where the potential 
V(R) has the minimum, I = (2 / pV") 1 / 4 is the amplitude 
of zero-point motion about this minimum, and e, given 
by 



e=\l 2 (p\V{R,)\ 



2 ) 



(B7) 
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is the corresponding energy in units of the oscillator fre- 
quency u> = 2 1 [it 1 . For the ground state we expect 



Comparing Eqs. (B5) and (BIO), we obtain 



< 1. 



(B8) 



The negative sign in the argument of D £ -i/2 in 
Eq. (B6) is chosen to obtain an exponentially decaying 
wave at large negative x, i.e., from the left turning point 
R_ and towards the origin. At large positive x, that is, 
at R — R + 3> I, both decaying and growing exponentials 
are present. At such x the wavefunction can be cast into 
the quasiclassical form 



V V2i 



(B9) 



which is equivalent to 



-^£±^2,/- C2 e- 2S ( b )- 2 . 
2Wire c_ V e 



(B12) 



For x at which the above calculation is valid we have 
S(b) ~ S - 1, where So = S(R = oo,x = 0). Thus, we 
arrive at 



»\V(Ro)\- ¥ 



fJ^e"-o 



(B13) 



Vl ip(R) 



-s(f>), 



i/>-(R) + c+e S{b >iP + (R) , (BIO) which leads to Eqs. (39)-(41) of Sec. II. 



see Eqs. (38), (B2), and (B6). This is our second con- 
nection formula except we still have to specify the pre- 
exponential factors c + and c_. In fact, only their ratio 
is important. With the help of the asymptotical expan- 
sion 28 for Dg, one finds it to be 46 



(Bll) 



Finally, a minor technical comment is in order. Since 
we have used the quasiclassical approximation, all coeffi- 
cients in Eq. (B13) have a relative accuracy O {e~ s ^). 
In particular, we expect that in place of V{Rq) we have a 
slightly more negative value, so that the ground-state en- 
ergy Eb never exceeds the oscillator ground-state energy 
V(Ro) + l/(2fil 2 ), as required by physical considerations. 
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